Displaying Geo-Spatial Data with R
Martin Schweinberger
2020-12-05
Introduction
This tutorial introduces geo-spatial data visualization in R. The entire R markdown document for this tutorial can be downloaded here.
This tutorial is based on R. If you have not installed R or are new to it, you will find an introduction to and more information how to use R here. For this tutorials, we need to install certain packages from an R library so that the scripts shown below are executed without errors. Before turning to the code below, please install the packages by running the code below this paragraph. If you have already installed the packages mentioned below, then you can skip ahead ignore this section. To install the necessary packages, simply run the following code - it may take some time (between 1 and 5 minutes to install all of the libraries so you do not need to worry if it takes some time).
# clean current workspace
rm(list=ls(all=T))
# set options
options(stringsAsFactors = F) # no automatic data transformation
options("scipen" = 100, "digits" = 4) # suppress math annotation
op <- options(gvis.plot.tag='chart') # set gViz options
# install libraries
install.packages(c("RgoogleMaps", "ggmap", "mapproj", "sf",
"dplyr", "OpenStreetMap", "devtools", "DT"))
# install package from github
devtools::install_github("dkahle/ggmap", ref = "tidyup")Depending on the maps that are used in the visualization, it may also be necessary to access other data bases. One very useful data base for maps is, of course, Google Maps. However, to access Google Maps materials, installation and setting up other pieces of software is necessary. How to get access to Google’s data is discussed below. In the following section, methods that do not require installation of software other than R.
1 Getting started with maps
The most basic way to display geospatial data is to simply download and display a map. In order to do that, we load the libraries necessary for extracting and plotting the map.
The package OpenStreetMap offers a range of maps with different features. To access the OpenStreetMap data base, it is necessary to install the package. Once the package is installed, we can simply extract the map and define the region we want to plot by defining the longitude and latitude of the upper left and lower right corner of the region we want to display. The argument minNumTiles defines the accuracy of the map, the higher the number of tiles, the higher the resolution. The type of map is defined by the type argument. The type argument defines from which server the map is extracted. Once we have extracted a map, we can plot it using the “plot” function.
# load library
library(OpenStreetMap)
# extract map
AustraliaMap <- openmap(c(-8,110),
c(-45,160),
# type = "osm",
# type = "esri",
type = "nps",
minNumTiles=7)
# plot map
plot(AustraliaMap)
In order to obtain different map types, we change the type argument. The following options are available for type:
# load package
library(DT)
opt <- c("osm", "osm-bw","maptoolkit-topo", "waze", "bing", "stamen-toner", "stamen-terrain", "stamen-watercolor", "osm-german", "osm-wanderreitkarte", "mapbox", "esri", "esri-topo", "nps", "apple-iphoto", "skobbler", "hillshade", "opencyclemap", "osm-transport", "osm-public-transport", "osm-bbike", "osm-bbike-german")
opt <- data.frame(opt)
datatable(opt, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")Unfortunately, not all options work. If they do not work, then an error message is shown telling us that the number of tiles is not supported.
We can zoom in or out by either changing the “zoom” or the “minNumTiles” arguments - in both cases, the higher the number, the more fine-grained the dispalyed map. Let’s check out some examples for maps of Queensland.
# extract map
queensland1 <- openmap(c(-8,135),
c(-30,160),
type = "osm",
minNumTiles=6)
queensland2 <- openmap(c(-8,135),
c(-30,160),
type = "esri",
minNumTiles=6)
# plot maps
par(mfrow = c(1, 2)) # display plots in 1 row/2 columns
plot(queensland1); plot(queensland2); par(mfrow = c(1, 1)) # restore original settings
The leaflet function from the leaflet package creates a Leaflet map widget using html-widgets. The widget can be rendered on HTML pages generated from R Markdown, Shiny, or other applications. The advantage in using this function lies in the fact that it offers very detailed maps which enable zooming in to very specific locations.
# load package
library(leaflet)
# load library
m <- leaflet() %>% setView(lng = 153.05, lat = -27.45, zoom = 12)
# display map
m %>% addTiles()Generating maps with rworldmap and base R
Another data base that is very useful when certain maps is the rworldmap package. The rworldmap package contains the shape files for countries but also more fine grained-shape files that display the states of selected countries. The most basic data, however, simply represents the shapes of the countries in the world.
Using the worldmap package has the advantage that one is not dependent on third parties and their servers but can operate within R without being denied access due to e.g. copy right issues or server maintenance.
# load library
library(rworldmap)
# get map
worldmap <- getMap(resolution = "coarse")
# plot world map
plot(worldmap, col = "lightgrey",
fill = T, border = "darkgray",
xlim = c(-180, 180), ylim = c(-90, 90),
bg = "aliceblue",
asp = 1, wrap=c(-180,180))
The basic map shown above can then be modified and enriched with color coding to convey geospatial data. The following shows how to customize the world map.
# load library
library(maps)
# plot maps
par(mfrow = c(1, 2)) # display plots in 1 row/3 columns
# show map with Latitude 200 as center
map('world', xlim = c(100, 300))
# add axes
map.axes()
# show filled map with Latitude 200 as center
ww2 <- map('world', wrap=c(0,360), plot=FALSE, fill=TRUE)
map(ww2, xlim = c(100, 300), fill=TRUE)
par(mfrow = c(1, 1)) # restore original settingsGenerating maps with rnaturalearth and ggplot2
We can also use the data provided by the rnaturalearth and the rnaturalearthdata and use ggplot function from the ggplot2 package as well as the sf package to create very nice visualizations of geospatial data. The advantage over using rworldmap and base R lies in the fact that the code is easier to interpret and the visualizations are more appealing.
# load packages
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
# load data
world <- ne_countries(scale = "medium", returnclass = "sf")
# gene world map
ggplot(data = world) +
geom_sf() +
labs( x = "Longitude", y = "Latitude") +
ggtitle("World map", subtitle = paste0("(", length(unique(world$admin)), " countries)"))
We can also easily zoom in on certain areas in the map using the coord_sf function and also prettify the map by adding some custom features such as a compass.
library(rgeos)
library(ggspatial)
# gene world map
ggplot(data = world) +
geom_sf() +
labs( x = "Longitude", y = "Latitude") +
coord_sf(xlim = c(100.00, 160.00), ylim = c(-45.00, -10.00), expand = FALSE) +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
theme_bw()
We will now customize these basic maps and add information to them.
2 Customizing Maps
Displaying basic maps is usually less interesting because, typically, we want to add different layers to a map. In order to add layers to a map, we need to have a shape file, i.e. a file which contains information about borders or locations that can then be displayed in different colors. In other words, we need to have a shape object to add information to the map.
# extract locations
world_points<- st_centroid(world)
# extract labels
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))
# generate annotated world map
ggplot(data = world) +
geom_sf(fill= "gray90") +
labs( x = "Longitude", y = "Latitude") +
coord_sf(xlim = c(100.00, 180.00), ylim = c(-45.00, -10.00), expand = FALSE) +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
coord_sf(xlim = c(100.00, 180.00), ylim = c(-45.00, -10.00)) +
theme(panel.grid.major = element_line(color = "gray60", linetype = "dashed", size = 0.25),
panel.background = element_rect(fill = "aliceblue")) +
geom_text(data= world_points,aes(x=X, y=Y, label=name),
color = "gray20", fontface = "italic", check_overlap = T, size = 3)
However, it is often the case that we want to add information that is not already available. Therefore, we load the airports data set which contains the longitude and latitude of airports across the world. We will then use this data to show the locations of airports across the globe.
# load data
airports <- read.delim("https://slcladal.github.io/data/airports.txt",
sep = "\t", header = T)
# inspect data
datatable(airports, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")To display the locations of airports on a map, we first plot the map and then add a layer of points to indicate the location of airports. In addition, the “plot” functions offers various arguments for customizing the display, e.g. by changing the backgroundcolor (bg), defining the color of borders (borders), defining the color of the shapes (fill and col).
# plot data on world map
plot(worldmap, xlim = c(-80, 160), ylim = c(-50, 100),
asp = 1, bg = "lightblue", col = "black", fill = T)
# add points
points(airports$Longitude, airports$Latitude,
col = "red", cex = .01)
It is, of course, also possible to highlight individual countries.
# create data frame with iso3 country codes and number of visits
countriesvisited <- data.frame(country = c("AUS", "JPN", "FIN",
"CZE", "POL", "AUT",
"USA", "GBR", "IRL",
"DEU", "DNK", "FRA",
"NDL", "BEL", "ESP",
"HRV", "SVN", "NOR",
"ITA", "HUN", "ROU",
"BGR", "GRC", "TUR",
"CHE", "ARE"),
visited = c(5, 1, 2, 1, 1, 3, 4, 4, 5, 11, 1, 1, 2, 2, 4, 4,
1, 1, 3, 1, 1, 2, 1, 1, 3, 2))
# inspect data
datatable(countriesvisited, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")# combine data frame with map
visitedMap <- joinCountryData2Map(countriesvisited,
joinCode = "ISO3",
nameJoinColumn = "country")## 25 codes from your data successfully matched countries in the map
## 1 codes from your data failed to match with a country code in the map
## 218 codes from the map weren't represented in your data# def. map parameters, e.g. def. colors
mapParams <- mapCountryData(visitedMap,
nameColumnToPlot="visited",
oceanCol = "azure2",
catMethod = "categorical",
missingCountryCol = gray(.8),
colourPalette = c("coral",
"coral2",
"coral3", "orangered",
"orangered3", "orangered4"),
addLegend = F,
mapTitle = "",
border = NA)
# add legend and display map
do.call(addMapLegendBoxes, c(mapParams,
x = 'bottom',
title = "No. of visits",
horiz = TRUE,
bg = "transparent",
bty = "n"))
It is, of course also possible to show only a part of the map by defining the x- and y-axes limits of the plot window.
# get map
newmap <- getMap(resolution = "low")
# plot map
plot(newmap, xlim = c(-20, 59), ylim = c(35, 71),
asp = 1, fill = T, border = "darkgray",
col = "wheat2", bg = "gray95")
# add points
points(airports$Longitude, airports$Latitude, col = "red", cex = .5, pch = 20)
This is of course also possible to show Australian airports.
# plot data on world map
plot(worldmap, xlim = c(110, 160), ylim = c(-45, -10),
asp = 1, bg = "azure2", border = "lightgrey", col = "wheat1",
fill = T, wrap=c(-180,180))
points(airports$Longitude, airports$Latitude,
col = "darkblue", cex = .5, pch = 20)
In addition to the location of airports, it is also possible to show how many flights arrive at an airport. As this information is not provided in the airport data, we load the routes data which contains information about the routes that airlines fly.
# read in routes data
routes <- read.delim("https://slcladal.github.io/data/routes.txt",
sep = "\t", header=T)
# inspect data
datatable(routes, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")To show the number of arrivals at an airport (only in terms of how mayn routes end at that airport), we extract the number of routes that end in each airport.
# load library
library(plyr)
# extract number of arrivals
arrivals <- ddply(routes, .(destinationAirportID), "nrow")
names(arrivals)[2] <- "flights"
# inspect data
datatable(arrivals, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")We can now merge the airports and the arrival data set to combine the information about the location with the information about the number of routes that end at each airport.
# create arrival table
airportA <- merge(airports, arrivals, by.x = "ID", by.y = "destinationAirportID")
# inspect data
datatable(airportA, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")This table allows us to plot not only the location of airports but also the number of routes that end there.
# get map
australia <- getMap(resolution = "low")
# plot data on world map
plot(australia, xlim = c(110, 160), ylim = c(-45, -10),
asp = 1, bg = "azure1", border = "darkgrey",
col = "wheat2", fill = T)
# add points
points(airportA$Longitude, airportA$Latitude,
# define colors as transparent
col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3),
# define size as number of flights div. by 50
cex = airportA$flights/50, pch = 20)
The same map can also be created using the ggplot2 package which offers a very high degree of flexibility and allows for easy customization.
# create a layer of borders
ggplot(airportA, aes(x=Longitude, y= Latitude)) +
borders("world", colour=NA, fill="wheat1") +
geom_point(color="blue", alpha = .2, size = airportA$flights/20) +
scale_x_continuous(name="Longitude", limits=c(110, 160)) +
scale_y_continuous(name="Latitude", limits=c(-45, -10)) +
theme(panel.background = element_rect(fill = "azure1", colour = "azure1")) +
geom_text(aes(x=Longitude, y= Latitude, label=City),
color = "gray20", check_overlap = T, size = 3)
In addition, it may be useful to overlay an area with shading to indicate the density of airports orflights.
# create map with density layer
ggplot(airportA, (aes(x = Longitude, y= Latitude))) +
borders("world", colour=NA, fill="antiquewhite") +
stat_density2d(aes(fill = ..level.., alpha = I(.2)),
size = 1, bins = 5, data = airportA,
geom = "polygon") +
geom_point(color="red", alpha = .2, size=airportA$flights/100) +
# define color of density polygons
scale_fill_gradient(low = "grey50", high = "grey20") +
theme(panel.background = element_rect(fill = "lightblue",
colour = "lightblue"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
# surpress legend
legend.position = "none",
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank()) +
geom_text(aes(x=Longitude, y= Latitude, label=Name),
color = "gray20", fontface = "italic", check_overlap = T, size = 2,
alpha = sqrt(airportA$flights)/10)
# load library
library(ggmap)
# define box
sbbox <- make_bbox(lon = c(115, 155), lat = c(-12.5, -42), f = .1)
# get map
ausbg = get_map(location=sbbox, zoom=4,
# source = "google",
source = "osm",
# source = "stamen",
color = "bw",
# color = "color",
# maptype="terrain")
# maptype="terrain-background")
maptype="satellite")
# maptype="hybrid")
# maptype="toner")
# maptype="hybrid")
# maptype="terrain-labels")
# maptype="roadmap")
# create map
ausbg = ggmap(ausbg)
# display map
ausbg +
stat_density2d(data = airportA, aes(x = Longitude, y= Latitude,
fill = ..level.., alpha = I(.2)),
size = 1, bins = 5, geom = "polygon") +
geom_point(data = airportA, mapping = aes(x=Longitude, y= Latitude),
color="gray20", alpha = .2, size=airportA$flights/20) +
# define color of density polygons
scale_fill_gradient(low = "grey50", high = "grey20") +
theme(axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.background = element_rect(fill = "aliceblue",
colour = "aliceblue"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
# surpress legend
legend.position = "none")
# define box
sbbox <- make_bbox(lon = c(152.8, 153.4), lat = c(-27.1, -27.7), f = .1)
# get map
brisbane = get_map(location=sbbox, zoom=10,
maptype="terrain")
# create map
brisbanemap = ggmap(brisbane)
# display map
brisbanemap +
geom_point(data = airportA, mapping = aes(x = Longitude, y = Latitude),
color = "red") +
geom_text(data = airportA,
mapping = aes(x = Longitude+0.1,
y = Latitude,
label = "Brisbane Airport"),
size = 2, color = "gray20",
fontface = "bold",
check_overlap = T) 
3 Color Coding Geospatial Information
We can also use color coding of countries to convey information about different features of countries such as their population size (or results of political elections, etc.).
We use the data provided in the rnaturalearthdata and rnaturalearth packages, which contain information about the population size of countries, to color code and thus visualize differences in population size by country.
# extract world data
world <- ne_countries(scale = "medium", returnclass = "sf")
# create cut off points in Population
world$pop_estC <- base::cut(world$pop_est,
breaks = c(0, 500000, 1000000, 10000000,
100000000, 200000000, 500000000,
10000000000),
labels = 1:7, right = F, ordered_result = T)
# load library
library(RColorBrewer)
# define colors
palette = colorRampPalette(brewer.pal(n=7, name='Oranges'))(7)
palette = c("white", palette)
# create map
worldpopmap <- ggplot() +
geom_sf(data = world, aes(fill = pop_estC)) +
scale_fill_manual(values = palette) +
# customize legend title
labs(fill = "Population Size") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
# surpress legend
legend.position = "none")
# display map
worldpopmap
We can also scale the colors and add a legend to assist with the interpretation of the map.
# start plot
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
# customize legend title
labs(fill = "Population Size")
It is relatively easy to combine color coding with points. However, it is also relatively easy to color countries if all countries have values. If this is not the case (as in the example below), we need to include an additional step so that countries that are not mentioned also receive a color coding. In this example, the data contain countries and cities that I have visited along with their latitude and longitude.
# load data
visited <- read.delim("https://slcladal.github.io/data/visited.txt",
sep = "\t", header = T)
# inspect data
datatable(visited, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")After loading the data, we determine how many cities I have visited in a given country and add this frequency to the data set.
# determine how often a country was visited
ncountry <- as.data.frame(table(visited$ISO3))
colnames(ncountry)[1] <- "ISO3"
# add frequency to visited
visited <- merge(visited, ncountry, by = "ISO3")
# inspect data
datatable(visited, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")The next part is tricky as we do not only need to determine the color for the countries I have visited but also the color for those that I have not visited. In order to do that, we load the data set that underlies the world map that will be displayed.
# load library
library(maptools)
# load world data for plotting
data(wrld_simpl)
# def. color (bias for contrast)
pal <- colorRampPalette(brewer.pal(6, 'Greens'),
bias = 10)(length(visited$Freq))
pal <- pal[with(visited, findInterval(Freq, sort(unique(Freq))))]
# define color for countries not in visited
countrycolor <- rep("white", length(wrld_simpl@data$NAME))
# define colors for countries in visited
countrycolor[match(visited$Country, wrld_simpl@data$NAME)] <- palAfter assigning colors to all countries, we can proceed by plotting the color coded map along with points for the locations of the cities I have visited.
# plot map
plot(wrld_simpl, ylim=c(-40, 85), xlim = c(-180, 180),
mar=c(0,0,0,0), bg="gray40", border = NA,
col=countrycolor)
# add points
points(visited$Longitude, visited$Latitude, col="red", pch=20, cex = .5)
A slightly more elegent way to achive the same goal is to plot the information using the ggplot2 package.
# use the contry name instead of 3-letter ISO as id
wrld_simpl@data$id <- wrld_simpl@data$NAME
wrld <- fortify(wrld_simpl, region="id")
# remove Antarctica from map
wrld <- subset(wrld, id != "Antarctica")
# change column names of visited
colnames(visited) <- ifelse(colnames(visited) == "Country",
"id",
colnames(visited))
# combine visited and wrld
worlddata <- join(wrld, visited, by = "id")
# change NA to 0
worlddata$Freq <- ifelse(is.na(worlddata$Freq) == T, 0, worlddata$Freq)
# convert freq into factor
worlddata$Freq <- as.factor(worlddata$Freq)
# start plotting
ggmapplot <- ggplot() +
geom_map(data=worlddata, map=worlddata,
aes(map_id=id, x=long, y=lat,
fill=Freq), color=NA, size=0.25) +
geom_point(data = visited,
aes(x = Longitude, y = Latitude),
col = "red", size = .75) +
scale_fill_manual(values=c("white", "gray80", "gray75",
"gray70", "gray65", "gray60", "gray55"),
name="No. cities visited") +
# coord_map("gilbert") + # spherical
coord_map() + # for normal Mercator projection
labs(x="", y="") +
theme(plot.background = element_rect(fill = "grey40", colour = NA),
panel.border = element_blank(),
panel.background =
element_rect(fill = "transparent", colour = NA),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = "right")
ggmapplot
4 Interactive Maps
An easy way to create interactice maps is to use the “googleVis” package. An example showing the population size by country provided by an in-built data set is shown.
# load library
library(googleVis)
# create motion chart object
Geo=gvisGeoMap(Population, locationvar="Country", numvar="Population",
options = list(height=350, daatMode='regions'))print(Geo, 'chart')A browser window should open and and if you accept and execute Adobe Flash, then an interactive map like the one above showing the population size by country. You should also be able to hoover over any country and see its population size.
You can also use the googleVis package to create customized maps. To customize an interactive map, we need to load a data set.
# load data
gapminderdata <- read.delim("https://slcladal.github.io/data/gapminder.txt",
header = T, sep = "\t")
# create map
Map <- data.frame(gapminderdata$country,
gapminderdata$year, gapminderdata$lifeExp)
# def. column names
names(Map) <- c("Country", "Year", "LifeExpectancy")
# inspect data
datatable(Map, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")We will only use the latest data (from 2007) and exclude all other data points and create the interactive map as we did above.
# remove data point not from 2007
Map <- Map[Map$Year == "2007",]
# create motion chart object
Geo=gvisGeoMap(Map, locationvar = "Country", numvar = "LifeExpectancy",
options = list(height = 350, dataMode = 'regions'))print(Geo, 'chart')Again, a browser window should open and and if you accept and execute Adobe Flash, then an interactive map like the one above showing the life expectancy by country. You should also be able to hoover over any country and see its available life expectancy.
5 Creating Maps with tmap and Shape Files
When creating more sophisticated maps, it is common to use shape files. Shape files define the edges of polygons and can be very complex depending on the number of edges of a polygon.
# load libraries
library(sf)
library(tmap)
# load shape file
ausstates <- sf::st_read("D:\\Uni\\UQ\\SLC\\LADAL\\SLCLADAL.github.io\\data\\shapes/AshmoreAndCartierIslands.shp", stringsAsFactors = F)## Reading layer `AshmoreAndCartierIslands' from data source `D:\Uni\UQ\SLC\LADAL\SLCLADAL.github.io\data\shapes\AshmoreAndCartierIslands.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 72.57811 ymin: -54.77708 xmax: 167.9966 ymax: -9.140693
## geographic CRS: WGS 84# plot australian admin shape file
tm_shape(ausstates) +
tm_fill() +
tm_borders()
# load library
library(rgdal)
# load shape files
ausstates <- readOGR(dsn = "D:\\Uni\\UQ\\SLC\\LADAL\\SLCLADAL.github.io\\data\\shapes/AshmoreAndCartierIslands.shp", stringsAsFactors = F)## OGR data source with driver: ESRI Shapefile
## Source: "D:\Uni\UQ\SLC\LADAL\SLCLADAL.github.io\data\shapes\AshmoreAndCartierIslands.shp", layer: "AshmoreAndCartierIslands"
## with 15 features
## It has 14 fieldsausstatesadmin <- readOGR(dsn = "D:\\Uni\\UQ\\SLC\\LADAL\\SLCLADAL.github.io\\data\\shapes/Australia_admin_6.shp", stringsAsFactors = F)## OGR data source with driver: ESRI Shapefile
## Source: "D:\Uni\UQ\SLC\LADAL\SLCLADAL.github.io\data\shapes\Australia_admin_6.shp", layer: "Australia_admin_6"
## with 252 features
## It has 68 fields
## Integer64 fields read as strings: z_order# plot australia and queensland
ozmap <- ggplot() +
geom_polygon(data = ausstates,
aes(x = long, y = lat, group = group),
colour = "black", fill = NA) +
geom_polygon(data = ausstatesadmin,
aes(x = long, y = lat, group = group),
colour = "red", fill = NA) +
theme_void()
# plot map
ozmap
# inspect shape file
summary(ausstates@data)## id country name enname
## Min. : 82636 NULL:AUS NULL:Ashmore and Cartier Islands NULL:Ashmore and Cartier Islands
## 1st Qu.:2316594 NULL:AUS NULL:Australian Capital Territory NULL:Australian Capital Territory
## Median :2316598 NULL:AUS NULL:Christmas Island NULL:Christmas Island
## Mean :2251865 NULL:AUS NULL:Cocos (Keeling) Islands NULL:Cocos (Keeling) Islands
## 3rd Qu.:2363491 NULL:AUS NULL:Coral Sea Islands Territory NULL:NA
## Max. :3225677 NULL:AUS NULL:Heard Island and McDonald Islands NULL:Heard Island and McDonald Islands
## NULL:AUS NULL:Jervis Bay Territory NULL:Jervis Bay Territory
## NULL:AUS NULL:New South Wales NULL:New South Wales
## NULL:AUS NULL:Norfolk Island NULL:Norfolk Island
## NULL:AUS NULL:Northern Territory NULL:Northern Territory
## NULL:AUS NULL:Queensland NULL:Queensland
## NULL:AUS NULL:South Australia NULL:South Australia
## NULL:AUS NULL:Tasmania NULL:Tasmania
## NULL:AUS NULL:Victoria NULL:Victoria
## NULL:AUS NULL:Western Australia NULL:Western Australia
## locname offname boundary adminlevel wikidata
## NULL:Ashmore and Cartier Islands NULL:NA NULL:administrative Min. :4 NULL:Q133888
## NULL:Australian Capital Territory NULL:NA NULL:administrative 1st Qu.:4 NULL:Q3258
## NULL:Christmas Island NULL:Territory of Christmas Island NULL:administrative Median :4 NULL:Q31063
## NULL:Cocos (Keeling) Islands NULL:Territory of the Cocos (Keeling) Islands NULL:administrative Mean :4 NULL:Q36004
## NULL:Coral Sea Islands Territory NULL:NA NULL:administrative 3rd Qu.:4 NULL:Q172216
## NULL:Heard Island and McDonald Islands NULL:NA NULL:administrative Max. :4 NULL:Q131198
## NULL:Jervis Bay Territory NULL:NA NULL:administrative NULL:Q15577
## NULL:New South Wales NULL:NA NULL:administrative NULL:Q3224
## NULL:Norfolk Island NULL:Territory of Norfolk Island NULL:administrative NULL:Q31057
## NULL:Northern Territory NULL:NA NULL:administrative NULL:Q3235
## NULL:Queensland NULL:NA NULL:administrative NULL:Q36074
## NULL:South Australia NULL:NA NULL:administrative NULL:Q35715
## NULL:Tasmania NULL:NA NULL:administrative NULL:Q34366
## NULL:Victoria NULL:NA NULL:administrative NULL:Q36687
## NULL:Western Australia NULL:NA NULL:administrative NULL:Q3206
## wikimedia timestamp
## NULL:en:Ashmore and Cartier Islands NULL:2018-09-11 02:23:50
## NULL:en:Australian Capital Territory NULL:2018-09-11 02:23:49
## NULL:en:Christmas Island NULL:2018-09-11 02:23:49
## NULL:en:Cocos (Keeling) Islands NULL:2018-09-11 02:23:50
## NULL:NA NULL:2018-09-11 02:23:50
## NULL:en:Heard Island and McDonald Islands NULL:2018-09-11 02:23:50
## NULL:en:Jervis Bay Territory NULL:2018-09-11 02:23:50
## NULL:en:New South Wales NULL:2018-09-11 02:23:48
## NULL:en:Norfolk Island NULL:2018-09-11 02:23:50
## NULL:en:Northern Territory NULL:2018-09-11 02:23:48
## NULL:en:Queensland NULL:2018-09-11 02:23:49
## NULL:en:South Australia NULL:2018-09-11 02:23:50
## NULL:en:Tasmania NULL:2018-09-11 02:23:48
## NULL:en:Victoria (Australia) NULL:2018-09-11 02:23:48
## NULL:en:Western Australia NULL:2018-09-11 02:23:48
## note rpath ISO3166_2
## NULL:NA NULL:2559345,80500,0 NULL:NA
## NULL:NA NULL:2354197,80500,0 NULL:AU-ACT
## NULL:NA NULL:2177207,80500,0 NULL:NA
## NULL:NA NULL:82636,80500,0 NULL:NA
## NULL:NA NULL:3225677,80500,0 NULL:NA
## NULL:NA NULL:2177227,80500,0 NULL:NA
## NULL:NA NULL:2357330,80500,0 NULL:NA
## NULL:NA NULL:2316593,80500,0 NULL:AU-NSW
## NULL:NA NULL:2574988,80500,0 NULL:NA
## NULL:NA NULL:2316594,80500,0 NULL:AU-NT
## NULL:NA NULL:2316595,80500,0 NULL:AU-QLD
## NULL:NA NULL:2316596,80500,0 NULL:AU-SA
## NULL:NA NULL:2369652,80500,0 NULL:AU-TAS
## NULL:border to be separated from riverbank (OSM community recommendation) NULL:2316741,80500,0 NULL:AU-VIC
## NULL:NA NULL:2316598,80500,0 NULL:AU-WAtable(ausstates@data$name)##
## Ashmore and Cartier Islands Australian Capital Territory Christmas Island Cocos (Keeling) Islands
## 1 1 1 1
## Coral Sea Islands Territory Heard Island and McDonald Islands Jervis Bay Territory New South Wales
## 1 1 1 1
## Norfolk Island Northern Territory Queensland South Australia
## 1 1 1 1
## Tasmania Victoria Western Australia
## 1 1 1ausstates_df <- broom::tidy(ausstates, region = "name")# extract names of states
cnames <- aggregate(cbind(long, lat) ~ id, data=ausstates_df, FUN=mean)
# define colors
# Define the number of colors you want
library(RColorBrewer)
nb.cols <- 15
newcolors <- colorRampPalette(brewer.pal(8, "Oranges"))(nb.cols)
# plot australian admin shape file
ausadminmap <- ggplot() +
# plot australian state shapes
geom_polygon(data = ausstates_df,
aes(x = long, y = lat, group = group, fill = id),
asp = 1, colour = NA) +
geom_text(data = cnames, aes(x = long, y = lat, label = id),
size = 2, color = "gray20", fontface = "bold",
check_overlap = T) +
# def. manual colors
scale_fill_manual(values = newcolors) +
# use empty theme
theme_void() +
scale_x_continuous(name="Longitude", limits=c(112, 155)) +
scale_y_continuous(name="Latitude", limits=c(-45, -10)) +
# def. background color
theme(panel.background = element_rect(fill = "slategray1",
colour = "slategray1",
size = 0.5, linetype = "solid"),
legend.position = "none")
# plot map
ausadminmap
# create data frame with longitude and latitude values
Latitude <- c(-25, -27.5, -25, -30, -30, -35, -25)
Longitude <- c(150, 140, 130, 135, 140, 147.5, 150)
mypolygon <- as.data.frame(cbind(Longitude,Latitude))
# inspect data
mypolygon## Longitude Latitude
## 1 150.0 -25.0
## 2 140.0 -27.5
## 3 130.0 -25.0
## 4 135.0 -30.0
## 5 140.0 -30.0
## 6 147.5 -35.0
## 7 150.0 -25.0# plot australian admin shape file
ausadminmap <- ggplot() +
# plot state shapes
geom_polygon(data = ausstates_df,
aes(x = long, y = lat, group = group),
asp = 1,
fill = "wheat1",
colour = "gray90") +
ggplot2::geom_text(data = cnames,
aes(x = long, y = lat, label = id),
size = 2, color = "gray20",
fontface = "bold",
check_overlap = T) +
# add customized polygon
geom_polygon(data=mypolygon,
aes(x=Longitude,y=Latitude),
alpha=0.5,colour="gray20",fill="red") +
ggplot2::annotate("text", x = mean(Longitude),y = mean(Latitude),
label = "My Polygon Area", colour="white", size=3) +
# def. manual colors
scale_fill_manual(values = newcolors) +
# add north arrow
ggsn::north(ausstates_df, scale = .1) +
# use empty theme
theme_void() +
scale_x_continuous(name="Longitude", limits=c(112, 155)) +
scale_y_continuous(name="Latitude", limits=c(-45, -10)) +
# def. background color
theme(panel.background = element_rect(fill = "slategray1",
colour = "slategray1",
size = 0.5, linetype = "solid"))
# plot map
ausadminmap
## Set options back to original options
options(op)Citation & Session Info
Schweinberger, Martin. 2020. Displaying Geo-Spatial Data with R. Brisbane: The University of Queensland. url: https://slcladal.github.io/maps.html (Version 2020.12.03).
@manual{schweinberger2020maps,
author = {Schweinberger, Martin},
title = {Displaying Geo-Spatial Data with R},
note = {https://slcladal.github.io/maps.html},
year = {2020},
organization = "The University of Queensland, Australia. School of Languages and Cultures},
address = {Brisbane},
edition = {2020/12/03}
}sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19041)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252 LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rgdal_1.5-18 tmap_3.2 googleVis_0.6.6 maptools_1.0-2 ggmap_3.0.0
## [6] ggspatial_1.1.4 rgeos_0.5-5 rnaturalearthdata_0.1.0 rnaturalearth_0.1.0 sf_0.9-6
## [11] maps_3.3.0 rworldmap_1.3-6 sp_1.4-4 leaflet_2.0.3 OpenStreetMap_0.3.4
## [16] hashr_0.1.3 stringdist_0.9.6.3 koRpus.lang.de_0.1-1 coop_0.6-2 hunspell_3.0
## [21] koRpus.lang.en_0.1-4 koRpus_0.13-3 sylly_0.1-6 textdata_0.4.1 magick_2.5.2
## [26] officer_0.3.15 gplots_3.1.0 FactoMineR_2.3 exact2x2_1.6.5 exactci_1.3-3
## [31] ssanv_1.1 ape_5.4-1 pvclust_2.2-0 NbClust_3.0 seriation_1.2-9
## [36] factoextra_1.0.7 cluster_2.1.0 Rmisc_1.5 plyr_1.8.6 psych_2.0.9
## [41] DescTools_0.99.38 pdftools_2.3.1 QuantPsyc_1.5 MASS_7.3-53 boot_1.3-25
## [46] car_3.0-10 carData_3.0-4 cfa_0.10-0 fGarch_3042.83.2 fBasics_3042.89.1
## [51] timeSeries_3062.100 timeDate_3043.102 Matrix_1.2-18 gutenbergr_0.2.0 sna_2.6
## [56] statnet.common_4.4.1 network_1.16.1 GGally_2.0.0 wordcloud_2.6 RColorBrewer_1.1-2
## [61] SnowballC_0.7.0 scales_1.1.1 likert_1.3.5 xtable_1.8-4 ggparty_1.0.0
## [66] RCurl_1.98-1.2 Boruta_7.0.0 Hmisc_4.4-1 Formula_1.2-4 survival_3.2-7
## [71] party_1.3-5 strucchange_1.5-2 sandwich_3.0-0 zoo_1.8-8 modeltools_0.2-23
## [76] randomForest_4.6-14 cowplot_1.1.0 Gmisc_1.11.0 htmlTable_2.1.0 Rcpp_1.0.5
## [81] flextable_0.5.11 caret_6.0-86 lattice_0.20-41 DT_0.16 partykit_1.2-10
## [86] mvtnorm_1.1-1 libcoin_1.0-6 vcd_1.4-8 knitr_1.30 kableExtra_1.3.1
## [91] gridExtra_2.3 tokenizers_0.2.1 tm_0.7-7 NLP_0.2-1 readxl_1.3.1
## [96] quanteda_2.1.2 tidytext_0.2.6 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
## [101] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.4 ggplot2_3.3.2
## [106] tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] class_7.3-17 foreach_1.5.1 lmtest_0.9-38 crayon_1.3.4 nlme_3.1-149
## [6] backports_1.1.10 reprex_0.3.0 rlang_0.4.8 rjson_0.2.20 glue_1.4.2
## [11] parallel_4.0.3 classInt_0.4-3 dotCall64_1.0-0 isoband_0.2.2 haven_2.3.1
## [16] tidyselect_1.1.0 usethis_1.6.3 rio_0.5.16 XML_3.99-0.5 ggpubr_0.4.0
## [21] magrittr_1.5 evaluate_0.14 RgoogleMaps_1.4.5.3 gdtools_0.2.2 cli_2.1.0
## [26] rstudioapi_0.11 rpart_4.1-15 fastmatch_1.1-0 sylly.en_0.1-3 leafem_0.1.3
## [31] fields_11.6 xfun_0.19 tmaptools_3.1 askpass_1.1 sylly.de_0.1-2
## [36] caTools_1.18.0 TSP_1.1-10 expm_0.999-5 ggrepel_0.8.2 png_0.1-7
## [41] reshape_0.8.8 ipred_0.9-9 withr_2.3.0 rle_0.9.2 bitops_1.0-6
## [46] slam_0.1-47 ranger_0.12.1 cellranger_1.1.0 e1071_1.7-4 pROC_1.16.2
## [51] coda_0.19-4 pillar_1.4.6 RcppParallel_5.0.2 multcomp_1.4-14 fs_1.5.0
## [56] scatterplot3d_0.3-41 raster_3.3-13 vctrs_0.3.4 ellipsis_0.3.1 generics_0.1.0
## [61] qpdf_1.1 lava_1.6.8.1 urltools_1.7.3 tools_4.0.3 foreign_0.8-80
## [66] munsell_0.5.0 compiler_4.0.3 abind_1.4-5 stars_0.4-3 rJava_0.9-13
## [71] forestplot_1.10 prodlim_2019.11.13 utf8_1.1.4 recipes_0.1.14 jsonlite_1.7.1
## [76] gld_2.6.2 lazyeval_0.2.2 latticeExtra_0.6-29 checkmate_2.0.0 rmarkdown_2.5
## [81] openxlsx_4.2.3 webshot_0.5.2 dichromat_2.0-0 rsconnect_0.8.16 yaml_2.2.1
## [86] systemfonts_0.3.2 htmltools_0.5.0 viridisLite_0.3.0 digest_0.6.27 assertthat_0.2.1
## [91] leaps_3.1 rappdirs_0.3.1 lwgeom_0.2-5 registry_0.5-1 spam_2.5-1
## [96] units_0.6-7 ggsn_0.5.0 mapproj_1.2.7 Exact_2.1 data.table_1.13.2
## [101] splines_4.0.3 labeling_0.4.2 broom_0.7.2 hms_0.5.3 modelr_0.1.8
## [106] colorspace_1.4-1 base64enc_0.1-3 mnormt_2.0.2 tmvnsim_1.0-2 inum_1.0-1
## [111] nnet_7.3-14 coin_1.3-1 fansi_0.4.1 ModelMetrics_1.2.2.2 R6_2.5.0
## [116] lifecycle_0.2.0 rootSolve_1.8.2.1 zip_2.1.1 curl_4.3 ggsignif_0.6.0
## [121] TH.data_1.0-10 iterators_1.0.13 gower_0.2.2 htmlwidgets_1.5.2 triebeard_0.3.0
## [126] crosstalk_1.1.0.1 janeaustenr_0.1.5 rvest_0.3.6 mgcv_1.8-33 flashClust_1.01-2
## [131] lmom_2.8 codetools_0.2-16 matrixStats_0.57.0 lubridate_1.7.9 leaflet.providers_1.9.0
## [136] gtools_3.8.2 dbplyr_2.0.0 gtable_0.3.0 DBI_1.1.0 httr_1.4.2
## [141] highr_0.8 KernSmooth_2.23-17 stringi_1.5.3 reshape2_1.4.4 farver_2.0.3
## [146] uuid_0.1-4 spatial_7.3-12 leafsync_0.1.0 xml2_1.3.2 stopwords_2.0
## [151] jpeg_0.1-8.1 pkgconfig_2.0.3 rstatix_0.6.0